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Presented in this paper is a technique that we propose for extracting the physical parameters 
of a rotating stellar core collapse from the observation of the associated gravitational wave signal 
from the collapse and core bounce. Data from interferometric gravitational wave detectors can 
be used to provide information on the mass of the progenitor model, precollapse rotation and the 
nuclear equation of state. We use waveform libraries provided by the latest numerical simulations 
of rotating stellar core collapse models in general relativity, and from them create an orthogonal 
set of eigenvectors using principal component analysis. Bayesian inference techniques are then 
used to reconstruct the associated gravitational wave signal that is assumed to be detected by an 
interferometric detector. Posterior probability distribution functions are derived for the amplitudes 
of the principal component analysis eigenvectors, and the pulse arrival time. We show how the 
reconstructed signal and the principal component analysis eigenvector amplitude estimates may 
provide information on the physical parameters associated with the core collapse event. 

PACS numbers: 02.70.Uu, 04.80.Nn, 05.45.Tp, 97.60.Bw 



I. INTRODUCTION 

The detection of gravitational radiation will likely 
come in the near future. LIGO P, [2, Q is at its ini- 
tial target sensitivity, and the detection of an event from 
an astrophysical source could come at any time. Around 
the globe a world-wide network of detectors is emerging; 
VIRGO in Italy GEO in Germany fe^ and TAMA 

in Japan p7j are operating alongside the LIGO detectors 
in the U.S. in the quest for gravitational wave detection. 
In a few years advanced LIGO and advanced Virgo will 
come on-line, with a sensitivity increase of 10 over initial 

LIGO 1,1, mini, increasing the prospective detection 

rate significantly. It will be a great day for physics when 
gravitational waves are finally directly detected, but it 
will also be the birth of a new way of observing the uni- 
verse and conducting astronomy and astrophysics. 

The observations of expected sources should be dra- 
matic: stellar collapse, pulsars, the inspiral of binary neu- 
tron stars followed by black hole formation, or even the 
stochastic background from the Big Bang. Gravitational 
wave astronomy will soon enter a regime where parameter 
estimation work will provide the means to make impor- 
tant astrophysical statements. Gravitational wave burst 
signals from rotating stellar core collapse and bounce, re- 
sulting in the formation of a proto-neutron star, are one 
of the more promising and potentially extremely valu- 
able sources for ground based interferometric detectors. 
Gravitational wave burst events are typically character- 
ized by their short duration (from a few milliseconds to 
about one second) and the absence of accurate theoret- 



ical predictions for their waveforms (as opposed to e.g. 
waveform templates for binary inspiral). If these events 
happen sufficiently close by, the ground based gravita- 
tional wave detectors will be able to observe them. Detec- 
tion searches for gravitational wave bursts typically use 
methods that can identify unmodeled but short duration 
events [ij. LIGO VIRGO 0, GEO [l|, 

and TAMA [ll[I3 have all recently conducted searches 
for gravitational wave bursts. 

The prediction of the exact burst signal to be expected 
from a rotating stellar core collapse event depends on the 
complex interplay of general relativity, nuclear and parti- 
cle physics. Furthermore, it is anticipated that the signal 
is produced by various emission mechanisms, first from 
the coherent motion of the collapsing and rebounding 
core during the proto-neutron star formation and then 
the ringing down of the nascent hot proto-neutron star 
(all in cases when the core rotates). This is then possibly 
followed by emission due to prompt convective motion 
behind the hydrodynamic shock in the central part of 
the star due to non-axisymmctric rotational instabilities 
for rapid rotation, as well as due to pulsations of the 
proto-neutron star, for instance, triggered by fall-back of 
matter onto it (see e.g. [18[ for a comprehensive review). 



Only recently, the full complexity of the prospective emis- 
sion mechanisms for gravitational waves in a stellar core 
collapse has become appreciated. Previously, research 
mainly concentrated on the signal contribution from the 
collapse, core bounce and ring-down phase, and there- 
fore this is the gravitational wave signal that is now by 
far most comprehensively and accurately investigated us- 
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ing numerical studies. Recently, a series of calculations 
of core collapse models with unprecedented physical re- 
alism has been performed [H, [lo, HH . These studies pro- 
vide strong evidence that the predictions of gravitational 
wave signals from the collapse and bounce phase arc now 
robust. The gravitational wave results from these inves- 
tigations have been made publicly available in the form 
of electronic waveform catalogs. If detected, such a gravi- 
tational wave signal can then ideally provide information 
on the mass of the progenitor model, the precollapse ro- 
tation, the density of the collapsed core, and the nuclear 
equation of state. 

However, even for the collapse and core bounce sig- 
nal (not to mention the possible signal contributions 
from other emission mechanisms in a stellar core col- 
lapse event) one cannot conduct a template based search 
as is done when looking for coalescing compact binary 
signals, since it would be computationally impossible to 
completely cover the signal parameter space, and there- 
fore various other methods have been developed to effi- 
ciently perform the search for gravitational wave burst 
signals in the detectors. Bayesian inferential methods 
provide a means to use data from interferometric grav- 
itational wave detectors in order to extract information 
on the physical parameters associated with an event ^22]. 
Markov chain Monte Carlo (MCMC) methods are a pow- 
erful computation technique for parameter estimation 
within this framework; they are especially useful in ap- 
plications where the number of parameters is large [23| |. 
Good descriptions of the positive aspects of a Bayesian 
analysis of scientific and astrophysical data are provided 
in [H, m, [H]. MCMC routines have been developed 
that will produce parameter estimates for gravitational 
wave signals from binary inspiral [H, |23, [H, IH, Is^, |3l| 
and pulsars (3^ . [ssl Is^ l . Accurate parameter estimation 
of observed gravitational radiation events is the neces- 
sary pathway to understanding important problems in 
astrophysics and cosmology. In this paper we present 
methods which are the starting point in the process of 
using MCMC methods to extract and estimate parame- 
ters associated with rotating stellar core collapse events 
from their gravitational wave signals emitted during the 
collapse, core bounce and ring-down phase. 

Unfortunately, due to the complexity of the event (in 
particular its intrinsic multi-dimensional nature), the 
computational time required to derive these signals in 
numerical simulations of rotating stellar core collapse is 
significant, and thus the waveform generation cannot be 
performed instantly while fitting a signal template to the 
data, and additional techniques to simplify the analysis 
are required. Instead of using waveforms corresponding 
to arbitrarily picked locations within parameter space, 
we reduce the complexity of the problem by simplifying 
the waveform space to the span of a small number of 
basis vectors. The basis vectors are derived from a rep- 
resentative catalog of numerically computed signal wave- 
forms through the use of principal component analysis 
(PC A) (ssj . The waveforms used in our present analysis 



are from the most recent, advanced, and comprehensive 
general relativistic study of rotating stellar core collapse 
by Dimmelmeier et al. [2l| and depend on the mass of 
the progenitor model, the precollapse rotation, and the 
nuclear equation of state. 

In this study we make use of a catalog that is smaller 
than eventually required in an extensive and accurate 
analysis of signals in real data, but with its non-trivial 
number of elements spanning a large portion of the in- 
teresting parameter space we expect it nevertheless to 
be sufficiently complex to probe and to demonstrate the 
method qualitatively. The MCMC algorithm essentially 
fits a superposition of derived basis vectors to the data, 
providing parameter estimates as well as associated un- 
certainties. We also demonstrate how PCA and MCMC 
then allow the reconstruction of the time series of a grav- 
itational wave burst signal, including confidence bands. 
We present initial results showing that these PCA basis 
vectors are actually physically meaningful, and that we 
can extract astrophysical parameters from the measured 
core collapse burst signal. The ultimate goal is to be able 
to extract information on the nuclear equation of state 
from the gravitational wave signal as it is observed by 
a network of interferometers, as well as other relevant 
physical parameters associated with the event, as far as 
this is possible from the collapse, bounce and ringdown 
signal alone. 

There have been other approaches to signal reconstruc- 
tion and parameter estimation with burst signals from 
stellar core collapse. The discipline essentially starts 
with the work of Giirsel and Tinto ^3^], who presented 
a method for reconstructing the time series for gravita- 
tional wave bursts. Rakhmanov [s^l developed a general 
scheme of Tikhonov regularization to be applied on a net- 
work of detectors in order to extract the signal. MCMC 
methods have also been applied to parameter estimation 
for burst events from cosmic string cusps [11] . Searle et 
al. 39*1 have recently worked out the Bayesian framework 
to the coherent detection and analysis of burst signals, 
including a wide range of special cases like white noise 
burst signals, but also signals constituting linear combi- 
nations of some set of basis vectors, like the ones used in 
this present study. Summerscales et al. [i^l proposed a 
maximum entropy technique in order to infer the time- 
dependent gravitational wave burst signals detected by 
a network of interferometric detectors; they then calcu- 
lated the correlation between their reconstructed wave- 
form and the entries in a catalog of rotating core collapse 
waveforms [4T'| in order to infer the physical parame- 
ters from the event. Our work is different in that we 
use waveforms from physical models (the table of rotat- 
ing core collapse signals) as our means of reconstructing 
the signal. PCA basis vectors are initially created from 
the table of waveforms, and the MCMC estimate of the 
amplitudes for the PCA basis vectors provides the re- 
constructed waveform. A linearly polarized gravitational 
wave signal can be reconstructed from the data of a single 
interferometer. The estimate of the physical parameters 



3 



comes from an association of the reconstructed waveform 
with the elements of the catalog. The results in this pa- 
per are derived only using data from a single detector 
(signals embedded in simulated noise matching that of 
LIGO at its target sensitivity). 

The organization of the paper is as follows. In Sec- 
tion|TT]we describe the methods by which LIGO, VIRGO, 
GEO, and TAMA are searching for burst events; the pa- 
rameter estimation technique we describe in this paper 
would be applied to event candidates from such a burst 
search. The method by which the catalog of stellar core 
collapse burst signals is established is described in Sec- 
tion mil A summary of the PGA method used in our 
study is presented in Section HVl while our analysis strat- 
egy is described in Section|Vl We illustrate how the anal- 
ysis method works when applied to simulated example 
data in Section IVIl and in Section IVIII we demonstrate 
that the PG basis vectors and the inferred corresponding 
coefficients actually carry physically meaningful informa- 
tion. In Section [Villi we provide a summary, and discuss 
the direction of future work. 



II. BURST SEARCH PIPELINES 

All collaborations associated with ground based inter- 
ferometric gravitational wave detectors have develo ped 
burst search pipelines in order to look for events [ij. 
Hi, m. Hi, nil. Although some gravitational wave burst 
sources are well modeled, the burst pipelines search for 
very general types of waveforms. The only assumption 
which is usually made concerns the duration of the sig- 
nal (less than a few hundred of ms) and the frequency 
bandwidth of the search, which can be as large as the 
detector bandwidth (from a few dozen of Hz up to the 
Nyquist frequency). Burst pipelines typically reduce the 
frequency bandwidth to a few kHz to focus on specific 
types of sources. 

Most of the gravitational wave burst pipelines look for 
an excess of energy in a time-frequency map using dif- 
ferent multi-resolution time-frequency transforms. An- 
derson et al. [4^ 1 initially considered the energy given by 
the Fourier transform in a frequency band. More recent 
methods [1^ make use of a wavelet decomposition of the 
data stream. In contrast, in [H, [4^, \^ a gravitational 
wave burst signal represented as a sine-Gaussian is im- 
plemeted. 

Usually the gravitational wave burst searches are per- 
formed using data from a network of detectors. Demand- 
ing that a gravitational wave burst event is seen simul- 
taneously in several detectors allows one to reduce the 
false alarm rate, which is rather high in a burst search 
due to the short duration of the signal. This is also the 
only way to disentangle a real gravitational wave sig- 
nal from transient noise events. There exist two kinds 
of network analysis: coincident or coherent filtering. In 
a coincidence analysis [H, \^ , each interferometer out- 
put is analyzed, providing a list of triggers. A coinci- 



dence in time and frequency is then required. A coherent 
analysis [sl, 113, IH, ill, \M, [U, [13] uses all interferom- 
eters' information by combining the input data streams 
or the filtered data streams into one single output which 
takes into account the detectors' antenna response func- 
tion assuming the source is at a given position in the 
sky. A coherent analysis more efficienctly suppresses non- 
stationarities that are expected to be incoherent in the 
different detectors. 

When using data from a network of detectors, burst 
pipelines can reconstruct the position in the sky of the 
source [H, [H, [H, [sHl . The accuracy of the source sky 
position depends on the time resolution of the pipeline. 
It is usually estimated for a set of different waveforms. 
Depending on the complexity of the signal waveform the 
time resolution can vary from a fraction of ms up to sev- 
eral ms (l2j . The frequency content (central frequency 
and bandwidth) of the event is also estimated by the 
burst pipelines. The frequency information can give some 
hints concerning the possible astrophysical source. For 
instance, a central frequency of a few hundred Hz would 
point towards the event of a core collapse in a massive 
star resulting in the formation of a proto-neutron star (2]| 
(in particular if the detection is accompanied by a neu- 
trino trigger following shortly after) while a higher fre- 
quency content at about 10 kHz could suggest that the 
event could be due to a neutron star collapsing to a black 
hole [5^. Besides, coherent pipelines can also extract 
from the data an estimation of the waveform without 
assuming any model for the triggers with sufficient am- 
plitude [Sy, |57[ . 

The work presented in this paper intends to extract 
even more information from the most significant events 
found by a gravitational wave burst pipeline assuming 
an astrophysical source model. The estimation of the 
parameters could also be used to reject a possible gravi- 
tational wave event candidate; the idea of distinguishing 
a real stellar core collapse signal from an instrumental 
glitch is briefly discussed in sections IVIII and IVIIII below. 
The starting point for our technique would be the list 
of candidate triggers produced by a burst all-sky search 
pipeline, or times provided by electromagnetic and/or 
neutrino observations (external triggers). In both cases, 
the MGMG will search for the event over some relatively 
small time span. The trigger times provide a relatively 
small number of data periods to be examined by our 
method. The MGMG would attempt to produce a re- 
construction of the signal based on a waveform catalog 
for stellar core collapse models. 



III. GRAVITATIONAL WAVE SIGNALS FROM 
ROTATING CORE COLLAPSE AND BOUNCE 

In stellar core collapse, the (possibly rotating) unstable 
iron core of a massive star at the end of its life contracts 
to supernuclear density on a dynamical time scale of the 
order of 100 ms. As the core material stiffens due to re- 
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pulsive nuclear forces, the collapse is halted abruptly, and 
the inner part of the core undergoes a rebound (the core 
bounce). After a period of ring-down oscillations, the hot 
proto-neutron star settles down, while the remainder of 
the star is possibly (mainly depending on the mass of the 
progenitor) blown off in a supernova explosion. 

If the precollapse stellar core is rotating, the dynam- 
ics of the evolution are reflected in the gravitational wave 
signal waveform with a slow rise during the core collapse, 
a large negative peak around bounce, and damped os- 
cillations in the ring-down phase. If rotational effects 
become important for rapidly rotating models, the core 
collapse can even be stopped by centrifugal forces alone 
at subnuclear density. However, in contrast to previous, 
less sophisticated studies of the stellar core collapse sce- 
nario (see e.g. [58| and references therein), Dimmelmeier 
et al. [20| have shown that the signal waveform remains 
qualitatively unaltered and is thus generic for a wide 
range of initial rotation strength. 

The gravitational radiation waveforms analyzed in this 
article are the burst signals from the most recent, ad- 
vanced, and comprehensive general relativistic study of 
rotating stellar core collapse models [21]. These simula- 
tions were performed with a computational code that uti- 
lizes accurate Riemann solvers to evolve the general rela- 
tivistic hydrodynamic equations and a nonlinear elliptic 
solver based on spectral methods for the fully cou- 
pled metric equations in the conformal flatness approx- 
imation of general relativity (60j . The precollapse iron 
core models were taken from recent stellar evolutionary 
calculations by Heger et al. [6l|, H^l , either with an intrin- 
sic or an artificially added rotation profile. These initial 
models were then evolved with a nonzero-temperature 
nuclear equation of state (EoS), either the one by Shen 



et al. IgsI. 

Swesty 



J (Shen EoS) or the one by Lattimer and 
6a | (LS EoS; with a bulk incompressibility 



K = 180 MeV), both in the implementation of Marek et 
al- [67j, including contributions from baryons, electrons, 
positrons, and photons. Deleptonization by electron cap- 
ture on nuclei and free protons during the collapse phase 
is realized as proposed and tested by Liebendorfer [68j . 

Of the 136 models investigated in the study by Dim- 
melmeier et al. [2l'|, 128 models have an analytic initial 
rotation profile. Their collapse behavior is determined 
by the following parameters: 

• The strength of rotation is specified by the prec- 
ollapse central angular velocity, which varies from 
ric,! = 0.45 to 13.31 rad (with the individual 
range depending on the differentiality of rotation). 
The infiuence of rotation strength on the collapse 
dynamics and on the burst signal via rotational flat- 
tening is very pronounced. 

For a wide range of slow to intermediate initial rota- 
tion strengths, the peak |/i|max of the gravitational 
wave amplitude is almost proportional to the ra- 
tio Tb/|VF|b of rotational energy to gravitational 
energy at bounce (which in turn increases approx- 



imately linear with f2c,i in this regime if all other 
model parameters are kept constant). 

For very rapid rotation, however, |/i|max reaches a 
maximum and declines again, when the centrifu- 
gal barrier slows down the contraction considerably, 
preventing the core from collapsing to high super- 
nuclear density. In such a case the frequency of the 
burst signal, which is practically constant for slow 
or moderate rotation, decreases significantly. Note 
that the trend to lower frequencies for very slowly 
rotating models shown in [2^, [2l| , which spoils the 
constancy of the signal frequency for such models, 
is due to a (possibly artificial) low-frequency con- 
tribution originating from post-bounce convection 
in the proto-neutron star, which superimposes the 
signal from the core bounce and proto-neutron star 
ring-down. 

• The difFerentiality of the precollapse rotation pro- 
file is set by a length scale with values ^ = 50 000, 
1 000, or 500 km, ranging from almost uniform to 
strongly differential rotation. For a fixed initial 
strength of rotation (as set by ric,i), a change in 
A alone has no strong effect on the gravitational 
wave signal, neither on the amplitude nor on the 
frequency. However, rapid rotation can only be 
achieved for comparably small values of A. 

• The set of precollapse cores encompasses models 
with a progenitor mass of Mpiog = 11.2, 15.0, 
20.0, or 40.0 Mq (masses at zero-age main se- 
quence). The choice of the progenitor mass has 
a direct impact on the mass of the inner core dur- 
ing the collapse and at bounce (and thus on the 
mass of the proto-neutron star). Already with- 
out rotation, the different progenitors produce an 
inner core at bounce with a mass that depends 
non-monotonously on the mass of the progenitor 
(with the inner core mass increasing in the order 
Mprog = 11.2, 20.0, 15.0, and 40.0 Mq). 

This variation is considerably amplified by rotation, 
which itself increases the mass on the inner core 
at bounce (approximately linear with ric,i at slow 
rotation and roughly quadratically at rapid rota- 
tion; see [2l|). Nevertheless, the reflection of that 
effect on the peak signal amplitude |/i|max is prac- 
tically negligible, as the effect of a high inner core 
mass (which causes a large quadrupole moment and 
thus a strong resulting gravitational wave signal) 
is canceled by strong centrifugal support, resulting 
in slower collapse dynamics and, consequently, a 
weaker gravitational wave signal. 

The main frequency of the gravitational wave burst 
signal is also not influenced significantly by the 
mass of the progenitor model, except that in mod- 
els with the least massive progenitor (Mprog — 
11.2 Mq) even strong initial rotation does not de- 
celerate the relatively small collapsing inner core 
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enough to make it bounce at sub nuclear densities. 
Hence the signal frequencies of such models re- 
main comparably high. We also point out that the 
approximate nature of the deleptonization scheme 
during the collapse phase up to core bounce em- 
ployed in these models could actually be responsi- 
ble for overemphasizing the variation in the mass of 
the inner core at bounce with respect to the mass 
of the progenitor. 

• The microphysical equation of state during the evo- 
lution is chosen to be either the Shen EoS or the 
LS EoS. In the subnuclear density regime, the two 
EoSs are rather similar, with the LS EoS being a bit 
softer on average. However, at supernuclear den- 
sity, the differences are more pronounced. Here the 
adiabatic index 7 of the LS EoS jumps to ~ 2.5, 
while in the Shen EoS 7 reaches values of ~ 3.0, 
making the nuclear material described by that EoS 
significantly stiffer. Consequently, the models with 
the Shen EoS consistently exhibit lower central den- 
sities at bounce and after ring-down. 

Still, with respect to the peak waveform amplitude 
l^lmax of the burst signal, this does not translate 
into an unequivocal trend due to a complicated in- 
terplay in the proto-neutron star between central 
compactness and density structure at intermediate 
radii. The peak frequency of the waveform spec- 
trum, on the other hand, is almost always lower 
if the Shen EoS is used, which is a direct conse- 
quence of the lower central densities in the bounce 
and post-bounce period due to the nuclear material 
stiffness from that EoS if all other parameters are 
identical and only the EoS is varied. 

In Fig. [1] we present a representative sample of wave- 
forms which illustrate the impact of a model's initial ro- 
tation state on the gravitational wave signal from core 
bounce, from slow and almost uniform initial rotation 
(model "s20alo05_shen"), moderately fast and modestly 
differential initial rotation (model "s20a2o09_shen"), 
to rapid and very differential initial rotation (model 
"s20a3ol5_shen" ). While the effect of varying the ini- 
tial rotation state is clearly reflected in the waveform for 
the models shown here, the influence of the progenitor 
model's mass and also the EoS for a fixed initial rotation 
state is much less apparent, and thus analogous pictorial 
waveform comparisons are not presented. 

We note that 8 progenitor models of [2l| are from a 
stellar evolutionary calculation that includes rotation in 
an approximate way. Hence, these models already have 
an initial angular velocity profile; no artificial rotation ac- 
cording to the simple analytic relation, like in the other 
models, is added prior to evolution. Consequently, as ini- 
tial rotation is not parametrized, these models are not of 
use for this study. However, the collapse dynamics and 
associated signal waveforms of these models are well rep- 
resented by models with artificially added precollapse ro- 
tation in terms of both signal amplitude and frequency. 
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FIG. 1: Sample waveforms of the core bounce signal for three 
models with varying initial rotation states while the mass of 
the progenitor model and EoS are fixed. Note the relatively 
small signal peak at the time of core bounce and the signif- 
icant late-time contribution from post-bounce convection for 
the slowly rotating model "s20alo05_shen" , and the overall 
lower signal frequency for the rotation-dominated and cen- 
trifugally bouncing model "s20a3ol5_shen" . The three sig- 
nals presented here adequately cover the waveform morphol- 
ogy of our signal catalog. 



and therefore it is justified to not separately consider 
their behavior here. The different influences of the vari- 
ous model parameters, which are summarized here only 
briefly, are discussed in detail in [2l| . The respective sig- 
nal data can be downloaded freely from an online wave- 
form template catalog 

We emphasize that the parameter selection for the 
models in [2l[ is fairly complete in that it accounts for all 
known relevant parameters which could have an impact 
on the burst signal from a core bounce (although, for in- 
stance in the case of neutrino effect, in an approximate 
way). Furthermore, the astrophysically meaningful range 
of the parameters has been reasonably exhausted (in the 
case of rotation, the expected strength of rotation for the 
bulk of stellar core collapse events is even at the lower end 
of the investigated range; see [t^). Only the selection of 
EoS is limited to two due to a lack of access to alterna- 
tive nonzero-temperature nuclear equations of state for 
stellar core collapse at the time when the study by Dim- 
melmeier et al. [21] was performed. However, with the 
Shen EoS and the LS EoS the two extremes of a rather 
stiff nuclear material and a somewhat soft one are prob- 
ably well covered. 

The range of frequencies and amplitudes of the signals 
is quite broad for the model sample of f2ll| , with the inte- 
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grated characteristic frequency fc spanning from about 
100 to 700 Hz (as opposed to the peak frequency /max 
of the spectrum which hes in a narrow interval around 
700 Hz; see Fig. 16 in 21]) and the integrated dimension- 
less characteristic amphtude ranging from 6 x 10^^^ 
to 7 X 10^^^ for e.g. initial LIGO and a distance to the 
source of 10 kpc (for the definition of these quantities, 
see the references in e.g. [2l|). 

Still, the main factor of influence causing this variation 
is the rotation of the core (as measured by r2c,i and A, 
and rather well resolved in the models of [21|), while the 
mass of the progenitor, precollapse differentiality of ro- 
tation, and microphysical EoS hardly affect both and 
he (see Fig. 17 in [21j). This partial degeneracy regard- 
ing the model parameters makes it very difficult to infer 
the unknown parameters from a detection of the gravita- 
tional wave burst signal emitted by a rotating stellar core 
collapse. In reverse, this means that parameters of little 
impact on the gravitational wave signal may be resolved 
only coarsely in a parameter study aiming at providing 
signal templates. 

For an event with at most moderate rotation, in prac- 
tice only the strength of rotation can be extracted from 
the waveform with confidence, provided the distance to 
the source and the orientation with respect to the rota- 
tion axis can also be determined. Only if a multitude 
of core collapse events can be detected via gravitational 
waves can systematic effects of e.g. the nuclear material 
stiffness (as described by the EoS) on the frequency also 
be analyzed with more certainty. 

While the independent parameters i}c,i, A, Mpiog, and 
EoS uniquely specify each collapse model, from the actual 
numerical simulation of the collapse a number of impor- 
tant quantities can be obtained which characterize the 
evolution of each individual model; these are typically 
presented as results for rotating stellar core collapse cal- 
culations. Among these are the maximum density /3niax,b 
at the time of core bounce, the ratio Tb/|W^|b of rota- 
tional energy to gravitational energy at bounce and the 
corresponding value Tpb/jW^lpb late after bounce. 

These quantities provide information about the col- 
lapse dynamics and, for instance, permit one to distin- 
guish between a core bounce that is mostly caused by 
the nuclear material stiffening at supernuclear densities 
to one that is dominated by centrifugal forces, which is 
again reflected by the waveform. Thus, in reverse the 
waveform encodes not only information about the inde- 
pendent model parameters, but also about these evolu- 
tion quantities. Consequently, in the analysis on the cor- 
respondence between PCs and physical parameters pre- 
sented in Section IVHl we not only consider the param- 
eters for the model setup but also the evolution of the 
quantities Pmax.b, Tb/|M^|b, and Tpb/jW^lpb- However, in 
contrast to these 'robust' parameters reflecting the global 
collapse dynamics we refrain from analyzing other quan- 
tities (like e.g. entropy at a specific off-center location or 
the time span of contraction) which are prone to depend 
more sensitively on the numerical evolution scheme or 



grid setup used for the model simulation. 



IV. SINGULAR VALUE DECOMPOSITION 

The signal waveforms used in the analysis reported 
here are originally generated at a higher sampling rate. 
Each waveform is subsequently resampled at a rate of 
16 384 Hz (the LIGO and GEO data samphng rate). The 
waveforms are then buffered with zeros so that they are of 
the same length. Finally, the waveforms are time shifted 
so that the first (negative) peak in each waveform (which 
occurs shortly after the time of core bounce) are aligned. 

For modeling purposes, the set of signal waveforms will 
be decomposed into an orthonormal basis. Following the 
method prescribed by [35| , we create a matrix H so that 
each column corresponds to a signal waveform from the 
catalog after subtracting the overall mean of the wave- 
forms. For m waveforms, each n samples long, H is a 
matrix with dimensions n x m. Using singular value de- 
composition [7l|, H is factorized such that 

H = USV', (1) 

where U and V are orthonormal n x r and ni x r ma- 
trices, respectively, and S is a diagonal r x r matrix 
containing the singular values of H in decreasing or- 
der, i.e. S = diag(si, Sr) with si > . . . > > 0, 
r = rank(H) < min(m, n). The columns of U are the 
eigenvectors of the empirical covariance matrix H H', and 
similarly, the columns of V are the eigenvectors of H'H. 
Additionally, the singular values in S are the square roots 
of the eigenvalues \i of either H H' or H'H, i.e. Si — \/Ai. 

The columns of U, i.e. iti, . . . , it^, form an orthonormal 
basis of the linear space spanned by the columns of H, 
i.e. the signal waveform space, and each signal waveform 
can now be uniquely represented as a linear combination 
of these eigenvectors. A measure of multivariate scatter 
about the mean is the trace of the empirical covariance 
matrix, tr(H H ), also called total variation^ which equals 
Y^i K- So the sum of the first k < r largest eigenvalues, 
-^ii measures how much of the total variability of 
the waveforms is explained by the first k eigenvectors. 
These are referred to as the first k "principal compo- 
nents" (PCs) and achieve an optimal dimension reduc- 
tion from the r-dimensional signal waveform space to a 
/c-dimensional subspace. 

Note that each waveform is n samples long and U has 
dimensions n x r. For the waveforms considered in this 
article, n is typically 1000 to 10 000 samples long, so com- 
puting the eigenvectors of H H in U can be computa- 
tionally expensive. On the other hand, the number m 
of waveforms in the catalog is of the order of 100. So, 
computing the eigenvectors of H'H in V is much less 
demanding, and they can be used to compute the eigen- 
vectors in U by first noting that 

H'H^;, = \,v,, (2) 
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FIG. 2: The top three principal components (PCs) derived 
from the catalog of waveforms described in Sec. IIIII 



where is the eigenvector of H'H corresponding to the 
ith largest eigenvalue A;. By pre- multiplying both sides 
with H, we obtain 



(3) 



Ui, i.e. Hi),; is the 



From this, we can see that Hvi 
eigenvector corresponding to the ith largest eigenvalue 
of H H' and thus equals the ith column of U. Fig. [2] 
illustrates the resulting first three PCs when applying 
this procedure to an actual waveform catalog. 



V. ANALYSIS STRATEGY 
A. Definitions 

The starting point of our analysis is a waveform cata- 
log, i.e. a catalog of time series (of equal lengths) describ- 
ing the gravitational wave signal of core collapse events 
corresponding to different input parameter settings. Let 
n denote the number of samples (discrete time points) 
in each waveform vector and m the number of waveform 
vectors in the catalog. As detailed in Section IIVI for 
these m time series we derive the first k < r eigenvec- 
tors, ail, ... , Xk, corresponding to the k largest eigenval- 
ues [sF]. Each of these eigenvectors again is of length n. 

The data to be analyzed (the noisy measurement) are 
given in the form of a time series vector y of length 
TV {N > n). This vector consists of additive non- 
white Gaussian noise with known (one-sided) power spec- 
tral density denoted by Si{f), superimposed by a core- 
bounce burst signal of length n located at an unknown 
instant T along the time axis. 

Our aim is to model the core-bounce burst signal ob- 
servation in terms of the basis of the k eigenvectors de- 



scribed above. To this end, we assume that the mean of 
the yi's is a linear combination PiXt^i + . . ■ + PkXi^k of the 
k eigenvectors, and zero before and after the burst signal. 
In matrix notation, this can be expressed in terms of the 
expectation value E[y] — X(7-)/3, where (3 = (/3i, . . . , Pk)' 
denotes the vector of regression coefficients, and the Nxk 
matrix -X^(t) has column vectors formed by the zero- 
padded k largest eigenvectors Xi, . . . ,Xk which are cycli- 
cally time-shifted by a lag T. 

The signal reconstruction is eventually accomplished 
based on the Fourier domain representations of data and 
signal; in the following we will be referring to the con- 
ventions explicated in Appendix |^ Let y denote the 
Fourier transformed data vector and a;^, i = 1, . . . , fc, de- 
note the discretely Fourier transformed eigenvectors after 
zero-padding each to length N . The real and imaginary 
parts of these form the columns of the N x k real-valued 
matrix X (neglecting the redundant elements due to her- 
mitian symmetry). 

One of the unknown parameters to be estimated is the 
signal's location T along the time axis, and in order to 
match data and signal, one needs to be able to shift both 
against each other in time. Let ^(t) be the matrix of 
Fourier domain basis vectors shifted in time so that these 
correspond to a particular signal arrival time T (with 
respect to some pivotal time point). Time-shifting of 
X by some lag T can be done directly in the frequency 
domain by multiplying the original Fourier transform by 
a factor of exp(— 27ri/T). 



B. Model 1: Basic linear regression model 

If the data were an exact linear combination of first k 
principal components plus measurement noise, then the 
following standard linear model would adequately model 
the situation: 



y = ^(T) /3 + e, 



(4) 



where e is the (Gaussian) noise vector with given (one- 
sided) spectral density Si{f). In the frequency domain 
this corresponds to 



y = ^(T) /3 + £, 



(5) 



where e now is the Fourier transformed noise vector. In 
Fourier domain, the real and imaginary components of 
the noise vector s then simply are independently zero- 
mean Gaussian distributed with variances proportional 
to the power spectral density Si{f ) [zl,!!^ (see also Ap- 
pendix [B| . 

What is known here are the data (measurement) y, the 
matrix X of basis vectors (derived from the waveform 
catalog through PC A), and the noise's power spectral 
density. The unknowns are the coefficients f3 and the 
time parameter T. The a priori information about these 
is expressed in the prior distribution P(/3,r), which is 
assumed to be uniform, i.e. any value is assumed equally 
likely. 
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C. Model 2: Random effects regression model 

In general, the basis vectors will not allow one to com- 
pletely reconstruct the original signal, there will always 
be a certain mismatch left unexplained when simplifying 
the problem to the reduced set of PCs as in ([5]) above. 
Neglecting the mismatch term would result in overcon- 
fidence in the reconstructed waveform and parameters. 
However, the simplified model introduced above will still 
be of interest, as it provides good approximations that 
are useful in the eventual implementation. 

Adding an extra mismatch error component m to the 
model, a random effect in statistical terminology, the 
model is now 

y = /3 + m + £. (6) 

In the following, we expect the original catalog to be 
sufficiently densely populated, so that the mismatch is 
primarily due to the effect of neglecting PCs, and not 
due to the signal being much unlike the ones in the cat- 
alog. All we assume to be known about the mismatch 
is that it contributes a certain fractional amount to the 
signal's power. That amount was expressed in [SSJ in 
terms of the match parameter fj,. Taking the difference 
between reconstructed and actual signal to be zero on 
average, and then expecting a certain (fractional) power 
for the mismatch, we use a Gaussian distribution with 
corresponding mean and variance for modelling the mis- 
match, 

PK|ai)=N(0,a2j ^ P(m,|a2j^N(0,fa^), (7) 

i.e. the mismatch m enters the model as an additional 
white noise component with a one-sided power spectral 
density of 2At(T^. The assumption of a Gaussian distri- 
bution is a simple and convenient choice here, and ac- 
cording to maximum entropy theory it also constitutes 
the most conservative possible choice [zl]. It also seems 
to perform well in practice, despite the fact that the ac- 
tual mismatch follows a rather heavy-tailed distribution 
with many near-zero and also a substantial number of 
extreme values. Since the mismatch is supposed to scale 
with the signal (i.e. the 'relative mismatch' is assumed 
constant), the mismatch's variance parameter cr^ is set 
such that 

a^«7^^||X/3f (8) 

depending on the PC coefficients /3 via the implied sum- 
of-squares (or power) of the PC contribution to the sig- 
nal, 

1 ^ N z 2 

^=l 3=1 

and the scaling factor 7^ = ^^t^ is set so that it corre- 
sponds to a particular match fi as in [ssj . The above rela- 
tionship between and 7^ results from assuming the PC 



and mismatch contributions to the signal s = X/3 + m 
to be (approximately) orthogonal: m _L X(3 (i.e. the 
mismatch is defined as what is not spanned by the set of 
PCs). 

The actual amount of mismatch is another unknown, 
and the (approximate) scaling of with the signal 
power as in Eq. ([5]) is ensured through the definition of 
the prior distribution, which is set up as 

P(/3,T,a^)=P(T)xP(£^. (10) 

=P(/3)xP(a^|/3) 

As in the previous Section fVBl the priors P(T) and P(/3) 
are set to be independent and uniform, but the condi- 
tional prior distribution P((t^,J^) is taken to be a scaled 
inverse distribution: 

P.,,.(a^|/3) = Inv-x^ (z., j'iW^m') (11) 

with degrees-of-freedom parameter v and scale param- 
eter (7^ ||X/3|p), so that the prior certainty in the 
scale of cr,^„|/3 is defined through ly. For example, a 
specification of 7^ = 0.1 and i/ ~ 3 implies for the 
prior that P(0.038 < crV(i||X/3||2) < 0.85) = 90% = 
P(0.73 < II < 0.98), and a conditional prior mean of 
n<^lM = ^l'i\\X(3r. Setting = yields the 
(improper) Jeffreys prior, which does not depend on its 
prior scale parameter [75|. The Inv-^^ distribution was 
chosen here because it constitutes the conjugate prior dis- 
tribution for this problem [t^, which makes it a 'natural' 
choice and makes the eventual implementation particu- 
larly simple, as will be seen in Section IV Dl below. The 
parameters 7^ and v may now be set so that the prior 
distribution reflects the reconstruction accuracy to be ex- 
pected from the given set of k basis vectors derived from 
the waveform catalog at hand. 

D. Monte Carlo integration 

Inference on waveforms and parameters usually re- 
quires integrating the parameters' posterior distribution, 
as one is interested in figures like posterior expectations, 
quantiles, or marginal distributions. These are here de- 
termined using stochastic (Monte Carlo) integration, i.e. 
by generating samples from the posterior probability dis- 
tribution and then approximating the desired integrals 
by sample statistics (means by averages, etc). The gen- 
eration of samples from the posterior distribution is done 
using Markov chain Monte Carlo (MCMC) methods, that 
is, by designing a Markov process whose stationary distri- 
bution is the posterior probability distribution of interest, 
and which may then be numerically simulated step-by- 
step to produce the desired posterior samples. The gener- 
ated samples then allow one to explore marginal or joint 
distributions of individual parameters, or of functions of 
the parameters as in the case of signal reconstruction, 
where the posterior distribution of -I- m is of in- 

terest, and T, (3 and m are random variables. 
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In the case of the basic hnear model (see Section IV Bl) . 
the set of unknown parameters consists of the vector 
of PC coefficients /3 and the time shift parameter T. 
An MCMC algorithm here may be implemented as a 
Gibbs sampler [75], since the conditional posterior dis- 
tribution P(/3|T, y) of coefficients /3 for a given time T 
is already known and easy to sample from (see Eq. 
below). In a Gibbs sampler, generating samples from 
the joint distribution is carried out by alternately sam- 
pling from the two conditional distributions P(/3|T, y) 
and P(T|/3,y) Samples from P(/3|r,y) can be gen- 
erated straight away, while a simple "Metropolis-within- 
Gibbs" step, i.e. a nested Metropolis sampler implemen- 
tation flB^ . may be utilized for sampling from from the 
(one-dimensional) distribution P{T\f3,y). 

For the basic linear regression model (see Section lVBI) . 
the conditional posterior distribution of the PC coeffi- 
cients (3 for a given time shift T is a multivariate Gaus- 
sian distribution: 



P(/3|T,y)=N(AT,ST), 



where 



(12) 



(13) 



At - {X-^T)D-'X^T))-'X^T^D-'y, (14) 

and D = diag(iT^(/j)) is the noise's covariance matrix, a 
diagonal matrix of the individual variances, as given in 
Eq. jBH) ill. 

Sampling from the posterior distribution works simi- 
larly for the extended model (see Section IV C|) , in addi- 
tion to (3 and T, the mismatch rfi and mismatch vari- 
ance (T^ need to be sampled from. Note that the condi- 
tional posterior distribution P(/3|T, rh, a-f^^,y) is not sim- 
ply multivariate Gaussian any more, since changes in (3 
lead to different prior density values (via its effect on 
the signal's sum-of-squares value). The expression in 
Eq. (|12p still is an excellent approximation and is useful 
for defining a proposal distribution within the Metropo- 
lis step of the algorithm. The conditional distribution of 
the mismatch vector rh is independent Gaussian: 

P(m|/3,r,af„,y)=N(/i,S), (15) 

where 



Aj.Re = R.e((y-^(T)/3)j) 



2 m 



2 m 



-2 



JV -2 N c ( f \ 
2 "rn 4At "^lUjV 



(16) 



(17) 



and analogously for the imaginary parts Im(mj). For the 
mismatch variance cr^ the conditional posterior is 

P(a^|y,/3,T,m) = Inv-x'(«;,s2), (18) 

where the degrees of freedom k and scale are 

(19) 
(20) 



1/7 



N, 
2 j_ 

JV 



(see also Appendix ICj). 

The MCMC sampler was eventually implemented as a 
Gibbs sampler, alternately sampling from the three con- 
ditional distributions of 



N 



1. (3,T\y,rh,al, 

2. m|y,/3,r,CT^, and 

3. al\y,l3,T,rh. 

The first step is done in a "Metropolis- within- Gibbs" 
step, using a symmetric proposal distribution for T, the 
approximated posterior of /3|T, . . . from Eq. for the 
corresponding f3 proposal, and accepting/rejecting as in a 
usual Metropolis-Hastings sampler based on correspond- 
ing posterior and proposal probability density values. 
Samples for the second and third step may be generated 
directly. 



VI. IMPLEMENTATION AND APPLICATION 
A. Setup 

In the following examples we are using the waveform 
catalog described in Section UTTl containing 128 gravita- 
tional radiation waveforms from rotating core collapse 
and bounce [2l|, [6^. The basis vectors to be used for 
signal reconstruction are generated through a PCA (see 
Section |IV| Utilizing the same code and general 

setup for calculating the models described in 0, [69j . 
we then compute three new rotating stellar core collapse 
models with input parameter values that did not appear 
in the original catalog (but lie within the range of the 
catalog parameter space) and their associated waveforms. 
The aim is then to reconstruct these waveforms using the 
PCA method. 

The simulated data used here is 1 s in length and sam- 
pled at 16 384 Hz, superimposed with (simulated) non- 
white Gaussian noise, and Tukey-windowed. The shape 
of the noise curve is here taken to correspond to LIGO 
at its initial sensitivity as stated in [t^; this is the def- 
inition that is also implemented in the LIGO scientific 
coUaboration's LSC algorithm fibrary (LAL) [t^. The 
noise's spectral density Si{f) is estimated by averaging 
over a thousand "empirical" periodograms of identically 
generated data (same size and resolution, same window- 
ing applied), more closely resembling a realistic case in 
which the noise spectrum would need to be estimated 
as well. In addition, this way convolution effects on the 
spectrum due to the finite data size and windowing are 
compensated. 

The prior distribution for the mismatch parameter ct^ 
is set by determining what match /x would be achievable 
for the given number k of PCs. This is done by project- 
ing each single waveform in the catalog onto the span 
of the PCs, and then (for any given number k of PCs) 
considering the distribution of achieved mismatch values 
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FIG. 3: Achievable match /i and scahng factor 7^ for given 
numbers k of PCs across the whole catalog used here. A 
particular number k of PC basis vectors yields a certain match 
fi for each of the 128 waveforms in the catalog. Shown here is 
the distribution of matches across the catalog for increasing 
values of k, as characterised by mean, median and some more 
quant iles. 



across the catalog. Fig. [3] illustrates the mean, median 
and some more quantiles of the mismatch distribution 
for increasing numbers of PCs (in analogy to Fig. 2 in 
Ref. fsHl)- We chose the prior's parameters (scale 7^ and 
degrees-of-freedom v) so that it matched the distribution 
of computed matches across the catalog. For the two ex- 
ample settings of fc = 10 and k — 20 considered in the 
following, this leads to maximum likelihood estimates of 
(72 = 13.3%, ly = 2.82) and (7^ = 5.70%, ly = 2.96). The 
resulting prior density for fc = 10 PCs is also shown in 
Fig.Sl 



B. Examples 

The first example illustrated here, labeled signal A in 
the following, has a signal-to-noise ratio (SNR) of p = 10, 
where the SNR is defined as 



\ 



j2 wWfjW 



It was generated assuming the Shen EoS, a mass Afprog = 
20 Mg for the progenitor, a precollapse central angular 
velocity fid — 5.48 rad s~^, a differentiality of the pre- 
collapse rotation profile given by A = 1000 km, and an 
effective distance of c?o = 5.17 kpc. The effective distance 
depends on the actual distance to the source, and also 
reflects the effect on the amplitude of the gravitational 
wave signal from other parameters (such as the geome- 
try of the detector with respect to the source), and is 
in general greater than the actual distance. The signal 
is embedded within simulated interferometer noise, and, 
using the model with A: = 10 basis vectors, the posterior 



distributions of the model parameters are derived using 
an MCMC implementation as described in the previous 
Section IV Dl The maximum achievable match for this 
example waveform (for fc = 10 and varying time shift T) 
is ^ = 0.97. 

The marginal posterior distributions of some individ- 
ual parameters are illustrated in Fig. |4l Note the timing 
accuracy, which has a standard error of 0.1 ms in this 
case; in the following examples (for the same SNR), these 
are also of the order of sub-millisecond. The parameters' 
joint posterior distribution provides the posterior proba- 
bility distribution of the waveform at any given time ti. 
Fig. [5] shows the posterior mean and a 90% confidence 
band in comparison with the originally injected wave- 
form. 

The example signal B illustrates how the signal recov- 
ery changes with differing numbers of PCs included in the 
model. This signal was generated assuming the LS EoS, 
a mass Mprog = 15.0 Mg for the progenitor, a precollapse 
central angular velocity fic.i — 2.825 rad s~^, a differen- 
tiality of the precollapse rotation profile oi A = 1000 km, 
and an effective distance of dc = 3.24 kpc. The signal's 
match is /i = 0.89 when using 10 PCs, and /i — 0.95 for 
20 PCs. The same data are analyzed twice, once using 
fc = 10 and then using fc = 20 PCs; the resulting recov- 
ered waveforms are shown in Fig. [Sj The enlarged model 
in principle allows for a better match of the signal, but on 
the other hand the larger number of parameters also de- 
creases the certainty in parameter estimates, so that the 
recovery does not necessarily improve. In this example, 
the posterior variances of the PC parameters common 
to both models (/3i, . . . ,/3io) as well as the time param- 
eter T increased for the fc = 20 case. In the resulting 
signal reconstruction (see Fig. [H]) the confidence band is 
wider, and the discrepancy between injected signal and 
posterior mean also increases. 

The reconstruction of signals at different SNRs can be 
seen in the following example signal C in Fig. [3 here the 
same signal is injected at different overall amplitudes and 
hence SNRs {p — 10 and p = 20). The injected signal 
was generated using the LS EoS, a mass Mprog = 40.0 Mg 
for the progenitor, a precollapse central angular veloc- 
ity ric,i = 9.596 rad s~^, a differentiality of the precol- 
lapse rotation profile given hy A — 500 km, and effective 
distances of dc = 8.32 and 4.16 kpc, respectively. The 
resulting signal waveform's match is p = 0.97. As ex- 
pected, a higher SNR yields a more accurate recovery, 
and the resulting posterior variances are correspondingly 
smaller. 



VII. CORRESPONDENCE BETWEEN 
PRINCIPAL COMPONENTS AND PHYSICAL 
PARAMETERS 

The analysis performed above not only allows one to 
reconstruct the waveform, but the posterior distribution 
of the PC coefficients actually also contains information 
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FIG. 4: Marginal posterior probability distributions of the parameters of example signal A with SNR p — 10, indicating the 
values and uncertainties of the model parameters as inferred from the data. Only the first four coefficients /9i to /34 corresponding 
to the top four principal components are shown here; there are fc = 10 coefficients in total. The dashed line in the top right 
plot shows the prior probability distribution in comparison to the (very similar) posterior. 
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FIG. 5: Reconstruction of the example signal A with SNR 
p = 10 using fc = 10 PCs. 



on the signal's physical parameters (ric,i, A, Mprog, and 
EoS) as well as other evolution parameters (like Pmax.b, 
^b/lW^lb, and Tpb/|PF|pb)- In the present study we are 
only able to choose between two specific EoSs, but as 
this work progresses the goal is to make estimates from a 
wider selection of EoSs. Considering the catalog of wave- 
forms alone, one can project each waveform onto the span 
of the PCs (as was done in the construction of Fig. [3]) and 
inspect the resulting fitted PC coefficients. Fig. [5] illus- 
trates the values of the first three coefficients {pi , P2 , and 
(3^) for all the signals in the catalog. Obviously the dis- 
tribution of coefficients in parameter space exhibits some 
structure. This structure can now be exploited, the idea 
being that the posterior means and variances of the coef- 
ficients Pi for a measured signal expose a pattern that is 
characteristic for the type of signal. For illustration, in 
the following we describe a simple approach involving the 
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FIG. 6: A repeated analysis and reconstruction of the same 
example signal B with SNR p = 10 using different numbers 
of PCs (fc = 10 and fc = 20). Note that a larger number of 
basis vectors in the model does not necessarily improve signal 
recovery. 



three example waveforms from the previous section. The 
coefficients' posterior distribution can be related back to 
the sets of coefficients for reconstructing every signal in 
the original catalog, in order to identify similarities. A 
simple ad hoc measure of similarity is a match that 
relates a set of best-matching coefficients Pi [i — 1, . . . , /c) 
associated with each catalog entry to the posterior dis- 
tribution of the coefficients V{Pi\y) of the inferred signal. 



12 




0.50 



0.52 
time (sec) 



0.54 



FIG. 7: Reconstruction of the same example signal C at dif- 
ferent overall amplitudes (and, with that, SNRs: p = 10 and 
p — 20), both times using fc = 10 PCs. 
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FIG. 8: The resulting first three PC coefficients /3i, (32, and 
133 when matching each of the 128 waveforms in the catalog 
against the set of PCs (no noise involved) . The waveforms in 
the catalog all correspond to a signal from a fixed distance 
(10 kpc). Note that the coefficients only occupy a very re- 
stricted region in parameter space. The diamonds show the 
posterior mean values of the PC coefficients for the three ex- 
amples (labeled A, B, and C) discussed in the text. For com- 
parison they are scaled down (by their known distance) to the 
same magnitude. 



Such a match may be defined as 
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considering only posterior means and variances of indi- 
vidual coefficients Pi, and allowing for a scaling factor c 



coefficient 

FIG. 9: Comparison of the distributions of all PC coefficients 
(3i to the best-matched sets of coefficients found in the origi- 
nal catalog. Posterior means and standard deviations ('error 
bars') are shown in black, the (scaled) sets of the top five best- 
matching coefficients are shown in different shades of gray. 



in the overall amplitude of the matched signal. Comput- 
ing the match for each signal from the catalog allows one 
to rank all signals and then determine those that fit the 
best. A mapping scheme set up this way combines ele- 
ments of nearest-neighbor and naive-Bayes classification 
techniques [t^ . 

For the initial example signal A shown above (see 
Figs. [4] and [5]), the posterior means and standard devia- 
tions ('error bars') are illustrated in a parallel coordinates 
plot [79] in Fig. ^ together with the five best-matching 
catalog entries. More details on the injected and the ten 
best-matching signals are given in Table H] Table |TT] lists 
the closest correspondents for the other example signals 
(examples B and C) in comparison with the injection 
values. In the case of signal A, the best matching wave- 
form from the catalog and the injected signal have very 
similar physical parameters. For signals B and C this 
simple method still provides a good estimation of the 
physical parameters, although admittedly not as accu- 
rate as case A; note that the reconstruction of the wave- 
form is quite accurate for all three signals. We think 
that through the use of techniques such as Procrustes ro- 
tation [iOl we will be able to provide a more direct link to 
the physical parameters. The posterior mean values for 
the first three PC coefficients (/3i to Ps) for all five exam- 
ples are also illustrated in Fig. [5] Reconstruction of the 
measured signals through the PC basis vectors not only 
allows one to capture the the waveforms' appearances, 
but also yields information on the underlying physical 
parameters. 

The information available from matching a signal can- 
didate against a set of basis vectors may not only facil- 
itate a classification within a set of possible astrophysi- 
cal sources (as illustrated in Fig. [5]), but may also help 
in distinguishing it from signals of different origin, like 
instrumental glitches, etc. Different types of potential 
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TABLE I: The injected example signal A and the corresponding top 10 best matching catalog entries. The column labeled 
with 'model' identifies the signal from the original catalog 69j. The following columns state achieved match (according 
to Eq. (|2H) ). the corresponding EoS, mass of the progenitor model (Mprog), precollapse angular velocity at the center (fici), 
precollapse differential rotation length scale (A), rotation rates initially (Ti/|W^|i), at the time of bounce (Tb/|W|b), and late 
after bounce (Tpb/IVKIpb), and the maximum density in the core at the time of bounce (pmax,b) ^3]. The corresponding 
best-matching effective distance de results from inverting the amplitude c (see Eq. (I2ip ') that yields the optimal match. 



rank 


model 




EoS 


-^-^prog 

[Mo] 


[rad s^^] 


A 

[km] 




b 

[%] 




Pmax,b 

[10" g cm^^*] 


dc 
[kpc] 


|W|T 1 


IW'Ipb 




example signal A (injected) 


Shen 


20.0 


5.48 


1000 


1.27 


11.9 


10.5 


2.62 


5.17 


1 


s20a2ol3_shen 


10.8 


Shen 


20.0 


6.45 


1000 


1.80 


14.8 


12.8 


2.42 


3.50 


2 


slla3ol2_shen 


12.8 


Shen 


11.2 


10.65 


500 


1.28 


14.9 


12.3 


2.46 


3.46 


3 


slla3ol3_shen 


13.7 


Shen 


11.2 


11.30 


500 


1.44 


16.1 


13.2 


2.36 


3.37 


4 


slla3ol3_ls 


13.8 


LS 


11.2 


11.30 


500 


1.44 


15.8 


13.2 


2.64 


3.10 


5 


s40a2o07_shen 


14.6 


Shen 


40.0 


3.40 


1000 


0.72 


11.8 


9.9 


2.58 


3.75 


6 


slla3ol2_ls 


15.4 


LS 


11.2 


10.65 


500 


1.28 


14.7 


12.3 


2.75 


3.12 


7 


sl5a2o09_shen 


15.4 


Shen 


15.0 


4.56 


1000 


1.09 


11.8 


10.3 


2.58 


3.42 


8 


s20a3o09_shen 


15.7 


Shen 


20.0 


8.99 


500 


0.90 


15.7 


12.5 


2.38 


3.58 


9 


s20a3o07_shen 


17.4 


Shen 


20.0 


5.95 


500 


0.50 


10.1 


8.0 


2.70 


3.24 


10 


s20a2ol3_ls 


20.6 


LS 


20.0 


6.45 


1000 


1.80 


14.4 


12.9 


2.75 


2.94 



TABLE 11: The injected signal example signals B and C, and the corresponding best matching catalog entries when varying k 
or p, respectively; see also Section [Vl Bl and Table HI 



PCs / SNR model 


EoS 


-^'-^prog 




A 


Ti 1 


Tb 






dc 








pmax,b 






[Mo] 


[rad s^"""] 


[km] 




[%] 


[10" g cm-=^] 


[kpc] 


example signal B (injected) 


LS 


15.0 


2.82 


1000 


0.41 


5.5 


5.2 


3.61 


3.24 


k = 10 slla3o07_shen 1.9 


Shen 


20.0 


5.95 


500 


0.40 


5.9 


5.1 


2.89 


3.41 


A: = 20 sl5alol3_ls 47.4 


LS 


15.0 


2.71 


50 000 


3.26 


6.1 


6.6 


3.56 


1.28 


example signal C (injected) 


LS 


40.0 


9.60 


500 


1.46 


22.1 


19.0 


0.73 


8.32 [4.16 


p = 10 s40a3ol3_ls 5.8 


LS 


40.0 


11.30 


500 


2.07 


23.4 


20.7 


0.28 


5.38 


p = 20 s40a2ol5_shen 19.6 


Shen 


40.0 


7.60 


1000 


3.62 


21.1 


20.3 


0.27 


2.32 



signals may also be partly reconstructible in terms of a 
linear combination of PCs, but a poor match, an miusual 
combination of coefficients, or a better match from an al- 
ternative set of basis vectors may then provide evidence 
in favor or against certain signal origins. 



VIII. CONCLUSIONS AND OUTLOOK 

In this paper we have described initial work on im- 
plementing a scheme by which the physical parameters 
associated with a rotating collapse and bounce in a stel- 
lar core collapse event can be estimated through the 
observation of gravitational waves by detectors such as 
LIGO, VIRGO, GEO or TAMA. We have presented a 
technique which allows one to reconstruct the detected 
signal through the use of numerically calculated gravita- 
tional radiation waveforms, principal component analy- 



sis, and Markov chain Monte Carlo (MCMC) methods; 
we have then compared the reconstructed signals to the 
table of numerically calculated waveforms, and inferred 
the physical parameters. 

The results displayed in this paper were achieved 
through the application of a number of different and dif- 
ficult analysis techniques. We are encouraged that in this 
initial study we have displayed that these methods can be 
combined in a way such that physical information about 
the supernova event can be extracted from the detected 
gravitational wave data. We have found that our method 
is quite successful in the reconstruction of the detected 
waveform, and this is displayed in the three examples. 
In the present analysis we have used the output from a 
single detector, and we expect the accuracy of the signal 
reconstruction will only get better as we expand the ap- 
plication of this work to the coherent analysis of multiple 
interferometers. We were able to make a fairly good asso- 
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ciation with the physical parameters; however, we think 
that we will be able to greatly improve on our ability 
to make parameter estimation estimates. In the present 
study we worked with a relatively small signal catalog, 
where each physical parameter took on only a handful of 
different values; for example, there were only two EoS to 
choose from. Our parameter estimation demonstration 
was pretty good within these constraints, but we expect 
the accuracy to increase with expanded signal catalogs. 
In addition, Procrustes rotation [80| is an example of one 
technique that we intend to apply in order to make bet- 
ter statistical estimates of the physical parameters based 
on the eigenstates provided by the PCA. 

Our eventual goal is to apply the method to the out- 
put from a network of detectors, which will provide an 
even better ability to discriminate signal and noise, re- 
construct the signal, plus allow us to estimate the source 
position on the sky. This would necessitate the use of 
catalogs for stellar core collapse waveforms which are 
both quantitatively more extensive and qualitatively im- 
proved, so that the parameter space of likely events can 
be completely covered and very finely resolved. For in- 
stance, a larger selection of available equation of state 
tables will extend the phase space in that respective di- 
rection, while a more accurate inclusion of neutrino ef- 
fects will both improve the quality of the burst signal 
from core bounce and as well yield a much more appro- 
priate characterization of the low-frequency signal from 
post-bounce convective bulk motion in the post-shock re- 
gion and in the proto-neutron star. In addition, other 
prospective emission mechanisms for gravitational waves 
in a core collapse event, like proto-neutron star pulsations 
or nonaxisymmetric rotational instabilities, could be also 
considered (for a comprehensive review on various such 
mechanisms, see (Tsj). 

The method presented in this paper will provide a 
better way to produce statistical and probabilistic state- 
ments about the physical parameters associated with a 
stellar core collapse event, and the ability to make state- 
ments about actual signals observed by interferometric 
detectors. The data from multiple interferometers will 
be coherently analyzed; this will allow the estimation 
of parameters associated with the location of the source 
on the sky, similar to what can be achieved with coher- 
ent MCMC analysis of binary inspiral gravitational wave 
data [1^, [sO]. As noted previously in our examples, the 
timing accuracy is better than 1 ms for a single interfer- 
ometer observation; a multiple interferometer coherent 
MCMC should provide very good sky localization esti- 
mation, and the use of pri ncip al components within the 
framework of Searle et al. [3^ also shows great promise. 
In addition, a coherent analysis will also allow us to bet- 
ter distinguish mismatch (m, common to all detectors) 
from noise (e, different between different detectors) and 
hence should improve the waveform reconstruction. We 
also plan to incorporate a flexible model for the detector 
noise spectrum, as described in [73| . A proper account- 
ing for the noise spectral densities of the interferometers 



could also be applied in producing the PCA eigenvec- 
tors and eigenvalues, i.e. these would not simply be the 
result of a least-squares fit, but rather a noise- weighted 
least-squares fit. The resulting basis vectors should then 
better reflect what is actually "visible" to an interfero- 
metric detector within its limited sensitivity band. 

An important part of our long-term research program 
will be to ensure that the catalog of waveforms truly 
spans and encompasses those from physically possible ro- 
tating stellar core collapses; it will be critical for us to 
consult closely with experts in the field J&, jM, HH to 
assure that the physical parameters for the waveforms we 
will use cover the range of natural possibilities. Typically 
when initially constructing the initial models one already 
has some idea about how many intermediate steps per 
free parameters are required. When the waveforms are 
generated as the output from the calculations (along with 
the other output data), then the variations of the wave- 
forms and the other data give some indication if the orig- 
inal choice of model construction and parameter division 
are sensible. While the signal catalog needs to stretch 
over the extreme limits of the values for the physical pa- 
rameters, a strength of our method is that the table need 
not be densely populated (as is the case with the tem- 
plate bank for coalescing binary signals); this the advan- 
tage provided by the use of the PCA. For example, when 
the PCA eigenvalues are derived, their variation provides 
guidance about the quality of our catalog. We can test 
whether the waveform table contains enough entries by 
successfully reconstructing additional (off-table) signals 
placed within the parameter space; this was essentially 
the method used for verifying the results presented in 
this article. 

The technique described in this paper could prove to 
be advantageous for signal reconstruction and parame- 
ter estimation when numerical techniques are required in 
order to produce the waveforms. For example, we can 
imagine our method as being useful for estimating pa- 
rameters for complicated binary inspiral signals; numer- 
ical relativity calculations are producing a wide array of 
complicated inspiral signals [8l| . We have already shown 
in Section IVIII that the fitted PC coefficients can be re- 
lated to the original physical parameters of the signals in 
the catalog. Future research will explore how the joint 
posterior distribution of the PC coefficients can be em- 
ployed for statistical inference on the physical parameters 
by linking the different parameter spaces via Procrustes 
rotation [80]. 

This method should also prove to be useful in distin- 
guishing a real stellar core collapse event from a com- 
mon noise glitch in the data. As displayed in Fig. [H 
there is a characteristic pattern to be found with PCA 
values associated with a real signal. A glitch produced 
by noise, when reconstructed via PCA, will likely fall 
outside the pattern formed by stellar core collapse sig- 
nals. Alternatively, one could set up an alternative glitch 
model, which, instead of using numerical simulations for 
the (PCA) basis vector generation, is trained on sets of 



15 



measured waveforms that are known to be instrumental 
artifacts. We intend to test this potential veto technique. 

In this paper we have reduced the table of 128 wave- 
forms to 10 or 20 of the most important PC eigenvec- 
tors. Instead of arbitrarily choosing the number of eigen- 
vectors to use, we plan to have our MCMC optimize 
the signal reconstruction through the use of a reversible 
jump MCMC, similar to the method described in [s^ . 
thereby treating the number or subset of PCs as another 
unknown. In this Bayesian approach the optimal num- 
ber of eigenvectors is determined by weighing the benefit 
of a good signal fit with a large number of eigenvectors 
against the Ockham's Razor penalization when the model 
gets overly complex. 
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APPENDIX A: DEFINITION OF DISCRETE 
FOURIER TRANSFORM 

The discrete Fourier transform is defined for a real- 
valued function h of time t, sampled at N discrete time 
points at a sampling rate of 

The transform maps from 

{h{t) e R : t = 0, At, 2At, ...,iN- 1)AJ (Al) 
to a function of frequency / 

{h{f) e C : / = 0, A/, 2A/, . . . , (A^ - 1)A/}, (A2) 
where A/ = and 

N-l 

Hf) - E Hj^t)exp{-27TijAJ). (A3) 

Since h is real-valued, the elements of h are symmetric 
in the sense that h{iAf) and h{{N — i)Af) are complex 



conjugates; this allows to restrict attention to the non- 
redundant elements indexed hy i — 0, . . . , N/2. 

APPENDIX B: FOURIER DOMAIN MODEL 

A Gaussian distribution of a (time-domain) random 
vector also implies gaussianity for the Fourier trans- 
formed vector. In particular, if n{t) is (zero mean, sta- 
tionary) Gaussian noise with one-sided power spectral 
density Si{f), then in the limit of large N and small 
At the real and imaginary components of the discrete 
Fourier transformed time series are independently zero- 
mean Gaussian distributed: 

P(Re(n(/,))) = N(0,a2^.), P(lm(n(/,))) = N(0,4), 
where the variance parameter is aj, — j^Si{fj) [t^-It^! 



APPENDIX C: JOINT AND CONDITIONAL 
POSTERIORS 

For the random effects model from Sec. IV CI the joint 
distribution of data and parameters may be factored out 
into: 



p(y,/3,r, CT^,jfi) 
= p(y|/3, T, al, rh) x p(m|/3, T, a^) 

xp(fT^|/3,T) xp(/3,T) (CI) 

ocexp(^-2g ) (C2) 

N 2 

xK)-"/%xp(-E^) (C3) 
x(^||X/3f)^(.f„)-^-^exp( ~"^^J^^^"' ),(C4) 

where the proportionality here refers to keeping the data 
and known parameters constant. In the above equation, 
the first term (|C2p is the 'usual' time-domain likelihood 
from Appendix iBl the term (jC3[) is the likelihood of the 
mismatch (random effect), and the term (IC4[) is a con- 
ditional prior of the mismatch variance parameter. By 
fixing the values for different subsets of parameters, one 
can determine conditional posterior distributions that are 
useful for the Gibbs sampling implementation (see Sec- 
tion |VDl). 
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